Reversible Gel Formation of Triblock Copolymers Studied by Molecular Dynamics 

Simulation 
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Molecular dynamics simulations have been employed to study the formation of a physical (ther- 
moreversible) gel by amphiphilic A-B-A triblock copolymers in aqueous solution. In order to mimic 
the structure of hydrogel- forming polypeptides employed in experiments [W.A. Petka et al, Science 
281, 389 (1998)], the endblocks of the polymer chains are modeled as hydrophobic rods repre- 
senting the alpha-helical part of the polypeptides whereas the central B-block is hydrophilic and 
semi-flexible. We have determined structural properties, such as the hydrophobic cluster-size dis- 
tribution function, the geometric percolation point and pair correlation functions, and related these 
to the dynamical properties of the system. Upon decrease of the temperature, a network structure 
is formed in which bundles of endblocks act as network junctions. Both at short and medium dis- 
tances an increased ordering is observed, as characterized by the pair correlation function. Micelle 
formation and the corresponding onset of geometric percolation induce a strong change in dynam- 
ical quantities, e.g., in the diffusion constant and the viscosity, and causes the system to deviate 
from the Stokes-Einstein relation. The dynamical properties show a temperature dependence that 
is strongly reminiscent of the behavior of glass-forming liquids. The appearance of a plateau in the 
stress autocorrelation function suggests that the system starts to exhibit a solid-like response to 
applied stress once the network structure has been formed, although the actual sol-gel transition 
occurs only at a considerably lower temperature. 

Keywords: molecular dynamics simulation, triblock copolymers, sol— gel transition, percolation transition 



I. INTRODUCTION 

Solutions of polymers with attractive groups (associ- 
ating polymers) exhibit a wide range of rheological prop- 
erties that can be controlled through variation of tem- 
perature and concentration. These materials enjoy ap- 
plications ranging from viscosity modifiers in food or oil 
recovery to adhesives and coatings (see Ref. p| and refer- 
ences therein). Under certain conditions, the attractive 
groups of the polymers associate to form a network and 
the system undergoes a sol-gel transition. The physical 
bonds between the attractive groups are reversible and, 
depending on their strength, can break and reform fre- 
quently on experimental time scales. The properties of 
these so-called weak or physical gels differ markedly from 
chemical gels in which the polymers are interconnected 
through covalent bonds. Compared to chemical gels, the 
current understanding of physical gelation is still lim- 
ited and even controversial Scenarios for thermore- 
versible gelation include the possibility of discontinuous 
gelation, in which the gelation is accompanied by sol-gel 
phase separation, and continuous gelation. The theoreti- 
cal treatments of Tanaka and Stockmayer 0, 0, 0, Q pre- 
dict that continuous gelation is a thermodynamic phase 
transition, whereas Semenov and Rubinstein arrive at the 
opposite conclusion 

Simulations can provide specific information that is 
not easily obtained otherwise. On the one hand, the mi- 
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croscopic structure of physical gels is difficult to deter- 
mine experimentally, whereas it can be directly accessed 
(within the limitations of the model used) by means of 
simulations. On the other hand, the numerical calcu- 
lations permit testing of theoretical hypotheses and ap- 
proximations. Two models are widely used to study asso- 
ciating polymers. In the first category of models, used in 
the above-mentioned theories, there are many association 
sites ("stickers") distributed along the polymer chain. 
This model was studied extensively in early (off-lattice) 
Monte Carlo simulations by Groot and Agterof 0, • 
Kumar and Panagiotopoulos have investigated the 
thermodynamic properties of a lattice-based version of 
this model by Monte Carlo simulations and did not find 
any indication that gelation is a thermodynamic phase 
transition. More recently, it has been observed that the 
dynamical properties of this model are similar to those of 
weak glass formers, in which the diffusion coefficient is de- 
scribed by an expression with a Vogel-Fulcher form |ll| . 

The second category of models consists of telechelic 
chains, in which the associative sites are located at both 
chain ends and typically represented by a single monomer 
(see Ref. 01 f° r a concise overview). At low concentra- 
tions, these chains have been predicted to form flower- 
like micelles ^3|. Simulations have indeed confirmed 
this Hi EH and found that the associative groups are lo- 
cated in the core of the micelles and the non-associative 
groups in the corona. At higher concentrations, micelles 
can be connected by "bridging" polymers for which both 
associative endgroups belong to different micelles, lead- 
ing to the formation of a micellar gel. The dynamics in 
such a system are governed by the hopping rate of the 
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associative groups between different micelies |l3f . Simu- 
lations have indeed found that the diffusion properties of 
such solutions can be described by an Arrhenius law 0] , 
as predicted by Tanaka and Edwards ^(|. If the poly- 
mer chains are less flexible, qualitative structural changes 
occur, as intra-chain pairing is suppressed and the forma- 
tion of flower-like micelles becomes energetically unfavor- 
able. This promotes the formation of a network structure 
at low polymer concentrations |T^ . However, the dynam- 
ical properties of such solutions of semiflcxible telechelic 
chains seem not to have been investigated. 

Recently, Petka et al. 0] have used genetic engineer- 
ing techniques to create artificial proteins consisting of 
a hydrophilic group flanked by two stiff hydrophobic 
blocks. This triblock copolymer was found to exhibit 
gelation in response to variation of pH or temperature. 
Its significance lies in the possibility to independently 
tune the strength of the endgroup attractions that are 
responsible for gelation and the solvent retention capa- 
bility of the chains, which is essential for the formation 
of a swollen gel. However, the actual structure of the hy- 
drogel, which is formed at low polymer concentrations, 
could only be conjectured. Motivated by the experimen- 
tal findings, we have employed molecular dynamics sim- 
ulations to investigate the dynamic and structural prop- 
erties of a solution of triblock copolymers that can be 
viewed as a greatly simplified, coarse-grained model of 
the artificial proteins. This model evidently does not 
capture all relevant properties of the experimental sys- 
tems, but rather should be viewed as a first attempt to 
determine the generic properties of a solution of triblock 
copolymers with two stiff endgroups. 



II. MODEL AND SIMULATIONAL DETAILS 

In order to study the gelation of triblock copolymers 
we employ molecular dynamics simulations, using the 
DL_POLY_2 code [f^. The polymers have an A-B-A 
structure, where the A-blocks are rigid hydrophobic rods 
and the B-block is hydrophilic and semiflexible. In our 
coarse-grained model, the solvent is modeled implicitly 
and each copolymer block is composed of spherical units 
("monomers") that represent an effective segment. The 
total length of each chain is set to 15 units, consisting 
of three A- monomers per hydrophobic block and nine B- 
monomers in the hydrophilic block. This choice is mostly 
based upon practical considerations. A minimum of three 
units is required to represent a rod-like endblock, whereas 
a longer central block would pose equilibration prob- 
lems, given the computationally accessible time scales. 
Monomers of type A interact via an attractive Lennard- 
Jones potential, 



pulsive, 



Uaa = 4eaa 



OA A 



(1) 



whereas the interactions between monomers of type B 
and the interactions between unlike pairs are purely re- 
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(2) 



We set £aa = £bb = £ab = £ and ctaa = obb = cab = c 
and cut off all interactions at 2.5<r. In order to express 
our results in reduced units we use e and a as units of 
energy and length, respectively. The reduced coupling, 
or inverse reduced temperature, J = e/k^T is varied be- 
tween 1 and 2 in the simulations. The semiflexible char- 
acter of the hydrophilic block is controlled by a harmonic 
angle-dependent potential, 



u e = hg(9-e ) 



(3) 



where kg = lOe/degree 2 and 8 = 175°. The value of 
#o was chosen in accordance with the model proposed 
in Ref. 01- I n combination with the large value for 
kg it causes the chains to adopt an extended structure 
(but without having the rod-like structure that would 
be obtained for 6 = 180°), and gelation is anticipated 
to occur at relatively low polymer concentrations. The 
a-helical structure in the artificial proteins |l8j is mim- 
icked by making the endblocks fully rigid. All monomer 
units within a chain are connected via a harmonic bond 
potential 



fcbond(?* 



r ) 



(4) 



where /cbond = 170e/<7 2 and ro = 1.30c 

The simulations are performed in the canonical (NVT) 
ensemble, in a cubic box of linear dimension L — 39a 
with periodic boundary conditions. The total number of 
chains equals N = 216, corresponding to a monomeric 
packing fraction of only 0.029, i.e., roughly twice the 
overlap threshold. The temperature is controlled by 
means of the Nose-Hoover thermostat |2jj. The equa- 
tions of motion are integrated using a "leap-frog" Ver- 
let scheme |2l[ . with a time step (in reduced units) 
At = 0.00287. In all runs, the system is first equili- 
brated for four million steps; for some low temperatures, 
even longer equilibration periods are used. Subsequently, 
40 million time steps are carried out for high tempera- 
tures and 200 million time steps for low temperatures. 
After the equilibration period, the configuration of the 
system is recorded every 1000 time steps for analysis of 
structural (e.g., chain conformations and percolation of 
the system) and dynamical (e.g., single-chain diffusion) 
properties. In addition, the energies and the stress ten- 
sors are calculated and recorded every 100 time steps for 
the calculation of the specific heat and the stress auto- 
correlation function. 

Despite the simplifications made in this coarse-grained 
model, the required simulation effort is still appreciable. 
The total amount of CPU time corresponds to approxi- 
mately 2.5 years on a single 2.0GHz Intel Xeon processor. 
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III. SIMULATION RESULTS 
A. Structural properties 

In order to characterize structural changes that take 
place in this system upon variation of the temperature, 
we employ an approach used in the study of micelle for- 
mation. For each configuration, bundles of endblocks are 
identified. An endblock is considered part of a bundle 
if its center monomer lies within a distance r c from the 
center monomer of an endblock that is already part of 
the bundle. Our results turn out to be insensitive to the 
precise value of r c 6 [a, 4a] and we have chosen r c = 2a. 
A configuration contains N(m) bundles of m endblocks 
and the bundle-size distribution is defined as the thermal 
average [l7j 



W{m) 



(N(m)) 



(5) 



Figure ^ shows W(m) for five values of the inverse tem- 
perature J, illustrating the formation of bundles of hy- 
drophobic blocks as the temperature is lowered. At high 
temperature (J = 1.00), the distribution function decays 
monotonically, with a single peak at to = 1 (isolated 
endblocks). This corresponds to a regular solution of 
chains that are not associated. For lower temperatures 
(J > 1.15), however, an additional "shoulder" appears in 
the distribution function, which develops into a peak that 
increases in height and shifts to larger bundle sizes if the 
temperature is further decreased. This signals the forma- 
tion of bundles in which large numbers of hydrophobic 
endblocks participate, resulting in important structural 
changes in the solution. The occurrence of an inflection 
point in W{m) has been taken as a criterion for the criti- 
cal micelle point ^tJ- Here, we associate the appearance 
of a shoulder in Fig. 2] with the onset of the bundling 
process. The corresponding characteristic inverse tem- 
perature J* lies between 1.15 and 1.20. 

A similar criterion was employed in Ref. |15| [but note 
that the quantity P(m) employed in this reference differs 
from lT"(m)], and a comparable temperature dependence 
was observed for the bundle-size distribution. There 
is, however, a marked difference in the morphology of 
the bundles. The flexible telechelic polymers studied in 
Ref. 0| form flower-like micelles in which the chains take 
the shape of a loop. Both endgroups of each chain lie in 
the core of the same micelle and the central block lies in 
the corona For more rigid chains, two effects occur. 
The stiff endblocks have a tendency to align inside the 
bundle, giving it the appearance of a microcrystalline do- 
main. In addition, the semiflexible character of the cen- 
tral block prevents the chain from adopting a loop-like 
conformation. This second effect is illustrated in Fig. |3 
in which we compare the average number of endblocks 
per bundle to the average number of endblocks in a bun- 
dle that belong to the same chain. As shown, single-chain 
loops are essentially absent at all temperatures, confirm- 
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FIG. 1: Bundle-size distribution [probability W(m) of en- 
countering a bundle containing m endblocks] for five represen- 
tative values of the inverse temperature J. The appearance 
of a secondary peak characterizes the formation of bundles in 
the solution. 
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FIG. 2: Comparison between the average number of end- 
blocks per bundle (closed squares) and the average number 
of endblocks per bundle that form a loop-like structure (open 
triangles). Whereas the average bundle size increases upon 
increasing coupling J (decreasing temperature), the average 
number of loops remains negligibly small. The rigidity of the 
chains thus prevents the formation of "flower-like" micelles in 
the solution. Error bars are smaller than the symbol size. 



ing the prohibitively large energy penalty incurred by 
ring formation. 

Instead, each chain takes an extended conformation, 
with both hydrophobic endgroups participating in dif- 
ferent bundles. Thus, even the formation of a continu- 
ous network becomes possible, in which the bundles of 
endblocks act as network junction points |17j . However, 
following the experimental observations in Ref. [lg, we 
have chosen a markedly lower concentration than in ear- 
lier simulation studies. Before investigating whether the 
formation of a network is nevertheless possible, we con- 
sider the energetic aspects of bundling. As illustrated 
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FIG. 3: The specific heat CV (in reduced units) as a function 
of inverse temperature J. The maximum is indicative of bun- 
dle formation in the system. The line serves as a guide to the 
eye. 




FIG. 4: Percolation probability of a solution of associative 
telechelic polymers at a monomeric packing fraction cj> = 
0.029, as a function of inverse temperature J. Since the per- 
colation probability equals unity for J > 1.15, we consider 
this to be a measure for the percolation threshold. 



in Fig. |21 the specific heat Cy exhibits a pronounced 
but relatively broad maximum around J = 1.23, corre- 
sponding to the creation of bundles of attractive end- 
blocks. The specific-heat maximum was found to occur 
at a temperature below the onset of micelle formation 
(at J = J*) in Ref. ^{|. Our data do not permit us to 
conclude this unambiguously. Indeed, for J > 1.20 the 
simulations become almost prohibitively expensive, ow- 
ing to the slow dynamic evolution of the system. Thus, 
the present system does not lend itself well to the appli- 
cation of finite-size scaling techniques for the determina- 
tion of the nature and precise location of the bundling 
transition. For example, in case of a continuous phase 
transition, the height of the specific-heat maximum will 
increase (up to corrections to scaling) as L a / V oc N a ^ 3 ^ . 
The exponent a/(iv) is typically rather small (e.g., 0.058 
for Ising-type criticality 22j), so that even doubling the 
number of chains would only increase the peak height by 
an amount comparable to the statistical accuracy of the 
data. 

In order to determine whether bundle formation indeed 
leads to the emergence of a connected network structure, 
we consider the percolation probability. Geometric per- 
colation of the polymer chains in the solution is a nec- 
essary condition for gelation. However, whereas chemi- 
cal gelation coincides with the occurrence of geometric 
percolation [2^ , physical gelation has been sug gest ed to 
take place only far below the percolation point [lj. We 
consider our polymer solution to be percolating if a con- 
nected path (composed of chains that bridge the bundles 
of endblocks) exists between any pair of opposite sides of 
the simulation cell. The percolation probability, which is 
defined as the probability that a configuration is perco- 
lating, is plotted as a function of inverse temperature in 
Fig- El The system always percolates for J > 1.15, i.e., 
near the characteristic inverse temperature J* for bundle 
formation. As the percolation probability certainly can 



exhibit strong finite-size effects, this determination must 
only be viewed as an estimate for the percolation thresh- 
old in the thermodynamic limit. The observation that 
J* and the percolation threshold coincide reaffirms our 
interpretation that the telechelic chains become intercon- 
nected through bundle formation and form an spanning 
network. This behavior appears to differ from what has 
been observed for the solution of flexible telechelic chains 
studied in Ref. |l5j | , which exhibits a comparable temper- 
ature dependence in the bundle-size distribution but is 
reported to exhibit geometric percolation at all tempera- 
tures. Figure [5] shows a typical configuration, obtained in 
a simulation performed at J = 1.30. A network of inter- 
connected hydrophobic junction points is indeed clearly 
discernable. 

The onset of percolation affects the single-chain con- 
formations as well. This is illustrated by means of the 
temperature dependence of the end-to-end distance i? , 
see Fig. For comparison, the figure also includes the 
end-to-end distance for an identical system in which the 
monomers in the endblocks experience a purely repulsive 
interaction [see Eq. (J5J]. Whereas R c increases for both 
systems as the temperature is lowered, the end-to-end 
distance increases more rapidly for the chains with at- 
tractive endblocks than for the purely repulsive chains. 
Because of the semiflexible character of the chains, the 
relative change in R D is only several percent, but never- 
theless the effect is clearly most pronounced for J > 1.15, 
i.e., near the percolation threshold J*. We ascribe it to 
the conformational changes induced by the network for- 
mation. Owing to the low polymer concentration, the 
bundles are relatively widely separated, forcing the con- 
necting chains to adopt an extended conformation. This 
observation is reinforced by considering the spatial cor- 
relations between endblocks. 

Indeed, the rod-like structure of the hydrophobic end- 



FIG. 5: Snapshot of a simulation at J = 1.30, for a system of 
linear size L = 78a (monomer packing fraction 0.029). The 
image represents approximately two-thirds of the simulation 
box. The hydrophobic endblocks are shown in red and the 
hydrophilic groups in cyan. The extended structure of indi- 
vidual chains as well as the bundling of hydrophobic blocks 
(cf. the peak in Fig. can clearly be seen. 



FIG. 6: Square of the end-to-end distance R B of the telechelic 
chains (closed triangles), as a function of inverse tempera- 
ture J. For comparison, this graph also shows the end-to-end 
distance for identical chains in which the endblocks do not 
possess an attractive interaction (open squares). For both 
chain types, R c increases upon increasing J (decreasing tem- 
perature), reflecting the decreasing flexibility of the center 
blocks. However, for J > J* ~ 1.15 the telechelic chains 
clearly exhibit a stronger tendency to adopt an extended 
structure, which is attributed to the formation of an inter- 
connected network. Error bars are smaller than the symbol 
size. 



blocks leads to an internal structure in the bundles that 
is absent in the models studied in Refs. 0, Il7| . Fig- 
ure {7\ shows the hydrophobe-hydrophobe radial distribu- 
tion function g(r) (calculated from their center-of-mass 
separation) at different values of J. As the temperature is 
decreased, two distinct features can be identified in this 
distribution function. The increasing maxima at short 
separations, which all lie at distances within the bundle 
size (cf. Fig. 0, correspond to intra-bundle alignment 
of endblocks. The emergence of this microcrystalline 
morphology can be understood from the fact that in an 
aligned bundle each endblock experiences a large num- 
ber of monomer-monomer interactions with surrounding 
endblocks. The cutoff distance employed in the Lennard- 
Jones potential is larger than the maximum distance 
between monomers on fully-aligned (close-packed) end- 
blocks, so that even a single pair of rods can have 9 pair 
interactions. As shown in Fig. El the average number of 
rods per bundle increases rapidly from approximately 2 
at J = 1.15 to almost 10 at J = 1.40, leading to tightly- 
bonded bundles. It is this bonding that makes the result- 
ing network resistant to external stress. A second feature 
arises in Fig.0at lower temperatures. As shown in the in- 
set, an additional peak appears at a position that roughly 
coincides with the calculated average bundle separation, 
which varies from 13 at J = 1.15 to 15 at J = 1.30. 
Thus, this peak characterizes the ordered arrangement of 
the bundles at low temperatures, and we conclude that 
the radial distribution function reflects the simultaneous 
emergence of both short-range and medium-range order 
upon cooling. 



B. Dynamical properties 

In order to determine whether the structural changes 
observed in the triblock copolymer solution indeed corre- 
spond to gelation, we consider the dynamical properties 
as a function of temperature. Evidently, bundle forma- 
tion and the formation of a percolating network are an- 
ticipated to have a strong influence on the diffusion prop- 
erties of the polymers. Figure |H1 shows the mean-square 
displacement of the center-of-mass of polymers at differ- 
ent values for J. At high temperatures, we observe the 
standard behavior in which the dynamics cross over from 
ballistic motion at short times to diffusive motion at long 
times. At low temperatures, an intermediate regime ap- 
pears where the dynamics are slowed down, indicative of 
the arrested dynamics resulting from network formation. 
Comparable observations were reported by Kumar and 
Douglas [llj] i n a Monte Carlo study of a lattice model of 
an associating polymer solution and by Bedrov el al. 
for micellar solutions, although it should be noted that in 
both studies the polymer concentration was considerably 
higher than in the current system (which has c/c* w 2) 
and that in the micellar system the change in dynamic be- 
havior was not associated with the formation of a network 
structure. The dynamic behavior seen in the low-density 
gel is similar to that found in glass- forming materials, 
but the underlying mechanism is different. The tempo- 
rary localization of the triblock copolymers is caused by 
the strong intra-bundle interactions experienced by the 
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FIG. 7: Hydrophobe-hydrophobe radial distribution function 
at different values of the inverse temperature J. At short 
distances, the existing peaks increase and new peaks emerge 
as J increases, an indication of the alignment of the rigid 
hydrophobic rods within a bundle toward a microcrystalline- 
like structure. At long distances, a new peak appears at the 
average bundle-to-bundle distance and grows higher at higher 
J, suggesting the structured arrangement of the bundles. 



endblocks, rather than by caging or jamming effects. We 
also note that the width of the ballistic regime depends 
on the simulation model. If water molecules are included 
explicitly or implicitly via a friction coefficient (Brown- 
ian dynamics), the ballistic regime may be rather narrow 
or even not observable at all. For all investigated tem- 
peratures, the polymers eventually diffuse and we extract 
the diffusion coefficient D from this long-time behavior, 
see Fig.^l At low temperatures, a second (narrow) diffu- 
sive regime arises for times between the ballistic and the 
intermediate regime. This corresponds to intra-domain 
motion of the endblocks. 

The diffusion coefficient exhibits an exponential depen- 
dence on inverse temperature over the entire temperature 
range that we investigated, but two regimes can be dis- 
cerned, separated near J* w 1.15. Since both percola- 
tion and micelle formation occur near this temperature, 
it is not possible to uniquely attribute the strong de- 
crease of D to either of these two phenomena. For the 
"sticker" model of Ref. [Tl| , in which association sites are 
distributed along the polymer chain, the diffusion coeffi- 
cient was found to be well described by a Vogel-Fulchcr 
law. In our system, the diffusive dynamics follow an Ar- 
rhenius law, D cx exp(— E/k-gT), where E is the effec- 
tive activation energy, similar to strong glass- formers (cf. 
Ref. |2J|) and micelle- forming telechelic polymers |15j . 
The Arrhenius-type behavior suggests that dynamic re- 
laxation, which takes place through the exchange of an 
endblock between two bundles (network junctions), is 
controlled by an energetic barrier. This barrier has a 
clear thermodynamic origin, namely the strong attrac- 
tion between endblocks. From Fig. E is estimated 
to be approximately lOe in the low-temperature regime, 
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FIG. 8: Mean-square displacement {r(t)) 2 of the single-chain 
center-of-mass, for inverse temperatures J = 1.00, 1.15, 1.20, 
1.24, 1.28, 1.32, 1.36, 1.40, 1.44 (only the lowest and highest 
value for J are labeled). In addition to ballistic motion at 
short times and diffusive motion at long times, a slow inter- 
mediate regime appears at low temperatures. 




FIG. 9: Diffusion coefficient D as a function of inverse tem- 
perature J on a log-linear scale. There are two regimes with a 
different exponential dependence on J, which are joined near 
the percolation point J* ~ 1.15. While the high-temperature 
data (J < J*) only permit an approximate fit, the diffusion 
coefficient in the low-temperature regime is clearly well de- 
scribed by an Arrhenius law, suggesting that activated pro- 
cesses control the relaxation of the system. Error bars are of 
the order of the symbol size or less. 



consistent with our earlier estimate of the number of in- 
teracting monomers in a pair of hydrophobic rods. Inter- 
estingly, the similarity between the dynamic properties 
observed in this system and those of the micellar system 
studied in Ref. [15| suggest that while the semiflexible 
character of the center blocks and the presence of rod- 
like hydrophobic endblocks change the structural proper- 
ties of the solution, these differences do not qualitatively 
affect the dynamic behavior. 
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FIG. 10: The viscosity r\ at different coupfings on a log-linear 
scale. Similar to D, r\ behaves differently in two regions di- 
vided at the percolation point J* = 1.15. In both regions r\ 
can be well described by an Arrhenius law. Scatter in the 
data at low temperatures is caused by uncertainties in the 
numerical integration ©. 



Since gelation will be accompanied by a dramatic in- 
crease in viscosity 77, we compute this quantity by inte- 
grating the stress autocorrelation function G(t) |2q. 



n 

Here, G(t) is defined as 
V 



G(t) dt 



G{t) = 



3k B T 



^2(<7 a p{to)Vat3{to +*)) 



(6) 



(7) 



a/3 



where V is the volume of the system and af3 assumes the 
values xy, yz, zx. The average (• • •) is taken over all time 
origins to. <Jaj3 is the stress tensor of the system pH 



N 
i=l 



1 N 

-Y 



fijft-Fijai 



(8) 



where m is the monomer mass, Vi a is the a-component 
of the velocity of atom i, rijp is the /3-component of 
the vector separating monomers i and j, and Fij a 
is the a-componcnt of the force exerted on monomer i 
by monomer j. The sum runs over all N monomers. For 
the calculation of G(t) we employ a fast Fourier trans- 
form |21|. which accelerates the calculation by several 
orders of magnitude compared to the direct calculation 
method. 

Figurc^|shows 77 as a function of J . Above the perco- 
lation threshold (J < 1-15), the viscosity increases grad- 
ually with decreasing temperature. However, in accor- 
dance with the behavior of the diffusion coefficient, 77 
starts to increase rapidly at the onset of percolation and 
micelle formation, and is described by an Arrhenius law. 




t 

FIG. 11: Normalized stress autocorrelation function G(t) for 
J = 1.22. The plateau value Go is small but distinctly 
nonzero, as shown in the inset. In addition, the plateau value 
does not exhibit a clear temperature dependence. 



Experimentally, gelation is characterized by the appear- 
ance of a plateau in the stress autocorrelation function 
G(t). In our simulations, we observe such a plateau for 
all temperatures below the percolation threshold. The 
plateau extends to longer times upon decreasing temper- 
ature, but eventually G(t) decays to zero. Since we use 
the "atomistic" (i.e., monomer-based) representation of 
the stress tensor a a p in Eq. J7J), the results exhibit rel- 
atively large fluctuations. Figure ^2 shows a representa- 
tive example. The plateau value Go is small, but clearly 
nonzero, as shown in the inset. Go is found to be only 
weakly dependent on temperature, and no clear trend 
can be identified, implying that the rapid increase in r\ 
(Fig. 1 1 U|l arises from an increase in relaxation time rather 
than from a variation in Go . The integration © partially 
suppresses the statistical fluctuations present in G(i), 
but the uncertainties in 77 still reflect the computational 
challenges, in particular at low temperatures. While the 
strong increase in viscosity follows unambiguously from 
Fig. ^3 we emphasize that the largest relaxation times 
(to be discussed below) are still much smaller than the ex- 
perimentally observed relaxation times for physical gels, 
which range from microseconds to seconds 1]. Thus, 
only for (computationally inaccessible) temperatures far 
below the percolation threshold would the system inves- 
tigated here undergo a sol-gel transition. 

Following Ref. |ll|. we employ the non-Gaussian pa- 
rameter OLi for the single-chain center-of-mass displace- 
ment m, 



a 2 



3 (rm 
5 (r{t) 2 ) 2 



1 



(9) 



to estimate a characteristic time. This parameter equals 
zero for both the ballistic and the diffusive regime. As 
shown in Fig. 1121 012 increasingly deviates from zero as 
the temperature is decreased, reflecting the heteroge- 
neous dynamics resulting from the hopping of endblocks 
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FIG. 12: The non-Gaussian parameter ati as a function of 
time at different inverse temperatures J. The deviation of Q2 
from zero is indicative of the heterogeneous dynamics caused 
by the hopping of endblocks between bundles. The peak po- 
sitions (squares) are used to define a characteristic time r. 




FIG. 13: The characteristic time r as extracted from Fig. 1121 
as a function of inverse temperature J on a log-linear scale. 
Two relaxation regimes can be discerned, separated by the 
percolation threshold J* ps 1.15. 
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between bundles. Also the maximum in Q!2, from which 
we extract a characteristic time r, shifts to larger times 
upon cooling. At short and long times, ai goes to zero, 
confirming the expected behavior for the ballistic and dif- 
fusive regimes, respectively. At intermediate times, how- 
ever, a2 deviates from zero and the deviation increases 
as the temperature is lowered. The positive deviation 
indicates the presence of anomalously fast chains, very 
similar to the heterogeneous dynamics observed in glass- 
forming liquids [23] • Figure displays r as a function 
of J, on a log-linear scale. We see that the values of 
r, which extend up to approximately 500to for the low- 
est temperatures, also follow an Arrhenius law, with a 
similar change in slope as observed for the diffusion co- 
efficient and the viscosity. The mean-square displace- 
ment corresponding to t (cf. Fig. [HJ is much smaller than 
the typical "hopping distance" or bundle separation (cf. 
Fig. UJ , comparable to what is observed for glassy sys- 
tems. Likewise, the observed values of r are much smaller 
than the relaxation times that follow from the extent of 
the plateau in G(t), which are approximately given by 
V/G (J). 

Finally, we demonstrate the similarity of our dilute 
polymer solution to a glass-forming liquid by considering 
how well the Stokes-Einstein relation is obeyed. For a 
Newtonian fluid, we expect the product of the diffusion 
coefficient and the viscosity to obey Dr] = k^T / '47ri?h, 
where i?h is the hydrodynamic radius. The product of D 
and rj is plotted as a function of J in Fig. [21 The be- 
havior of the end-to-end distance (Fig. [fjj) suggests that 
Rh will change only weakly with temperature, and even 
may increase as the temperature is lowered. In combi- 
nation with the linear temperature dependence of the 
numerator of the Stokes-Einstein relation (i.e., inverse 
dependence on J), the steep increase in Drj for J > J* 
indicates a clear breakdown of this relation once micel- 



0.1 1 — 1 — 1 — 1 — 1 — ' — • — ' — • — 1 — 1 — 1 — 1 — 1 — 1 — 1 
1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 

J 

FIG. 14: The product Dr\ as a function of J. Below J* = 1.15, 
Dr) varies weakly, consistent with the prediction of the Stokes- 
Einstein equation. Above J* , an order of magnitude increase 
of Dn within a narrow temperature windows clearly deviates 
from the Stokes-Einstein equation. 

lization and network formation set in. The diffusion co- 
efficient is larger than would be predicted on the basis of 
the viscosity shown in Fig. 1101 The precise microscopic 
origin of this behavior, however, remains to be deter- 
mined. As shown in the inset of Fig. ^] this breakdown 
does not result from a temperature dependence in the 
modulus Go- 



IV. CONCLUSION 

We have studied the structural and dynamical prop- 
erties of a solution of associative A-B-A triblock copoly- 
mers. The semiflexible character of the center (B) blocks, 
in combination with the attraction between endgroups, 
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allows these polymers to form a gel at remarkably low 
concentrations. The molecular dynamics simulations pre- 
sented here form a natural extension of earlier work on 
scminexible chains |Tt| . which however only addressed 
structural properties. In addition, we observe dynamic 
effects that bear close resemblance to those reported in 
Ref. f° r a lattice-based model studied by Monte Carlo 
simulations and to those reported in Ref. jl5| for micelle 
formation. However, the solutions in either of these stud- 
ies had a significantly higher polymer concentration. Fur- 
thermore, the micelle solution was found to exhibit geo- 
metric percolation at all temperatures and the dynamic 
changes were linked to the thermodynamic micelle tran- 
sition. We observe that all dynamic changes are corre- 
lated with micellization and the simultaneous emergence 
of a percolating network of polymer chains, in which bun- 
dles of rigid endblocks act as network junctions. Upon a 
further decrease in temperature, the hydrophobic blocks 
tend to align within a bundle, forming a microcrystallinc- 
like structure. The resulting strong binding of the chains 
is responsible for the mechanical stability of the gel-like 
network. On a larger scale, the bundles distribute more 
regularly at lower temperature, as indicated by the ap- 
pearance of peak at the average bundle separation in the 
hydrophobe-hydrophobe radial distribution function. 

The change in the dynamical behavior of the solu- 
tion that occurs upon micellization and network forma- 
tion is reflected in the diffusion constant, the viscosity 
and the maximum in the non-Gaussian parameter. The 
temperature dependence of all these properties changes 
near the percolation threshold and is well described by 
an Arrhenius law, similar to what is observed for strong 
glass formers. The activation barrier in our system has 



a clear thermodynamic origin, namely the strong attrac- 
tion between endblocks that are part of the same mi- 
crodomain or "bundle." Thus, there are similarities with 
the diffusion of diblock copolymers in a lamellar phase, 
which also exhibits an exponential decay with temper- 
ature [2^, |2^, |3(j ■ Recently, this type of dynamics has 
attracted attention in the context of slow dynamics in 
systems with frustration-limited domains (see Ref. |3l| 
and references therein). 

Finally, we note that, while we observe a finite plateau 
in the stress autocorrelation function, as would be ex- 
pected for a gel-forming material, the dynamics in our 
systems are still faster than in actual experimental gels. 
Thus, the sol-gel transition only occurs at a temperature 
far below the percolation threshold. 
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